function plt_D0_orbital(orb_center,orb_orientationx,orb_orientationz,wfsign)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here


scan1=linspace(0,pi,101);
scan2=linspace(0,2*pi,121);

[mesh1,mesh2]=meshgrid(scan1,scan2);

wf_dxy=(sin(mesh1).^2).*sin(2*mesh2);
wf_dx2y2=(sin(mesh1).^2).*cos(2*mesh2);

wf_dxz=2*(sin(mesh1).*cos(mesh1)).*cos(mesh2);
wf_dyz=2*(sin(mesh1).*cos(mesh1)).*sin(mesh2);

wf_dz2=(1/sqrt(3))*(3*cos(mesh1).^2-1);



% 
% 
% wf_px=sin(mesh1).*cos(mesh2);
% wf_py=sin(mesh1).*sin(mesh2);
% wf_pz=cos(mesh1);

shift_center=orb_center;
%wf_porb=orb_orientation(1)*wf_px+orb_orientation(2)*wf_py+orb_orientation(3)*wf_pz;

wf_porb=wf_dz2*wfsign;


wf_porb_pos=wf_porb;
wf_porb_pos(wf_porb_pos<0)=0;
wf_porb_pos=abs(wf_porb_pos).^2;

wf_porb_plt=wf_porb_pos;
pltx0=sin(mesh1).*cos(mesh2).*wf_porb_plt;
plty0=sin(mesh1).*sin(mesh2).*wf_porb_plt;
pltz0=cos(mesh1).*wf_porb_plt;


newx=orb_orientationx/norm(orb_orientationx);
newz=orb_orientationz/norm(orb_orientationz);
newy=cross(newz,newx);
pltx=pltx0*newx(1)+plty0*newy(1)+pltz0*newz(1)+shift_center(1);
plty=pltx0*newx(2)+plty0*newy(2)+pltz0*newz(2)+shift_center(2);
pltz=pltx0*newx(3)+plty0*newy(3)+pltz0*newz(3)+shift_center(3);


mapC=zeros(size(plty,1),size(plty,2),3);
mapC(:,:,1)=1;
%mapC(:,:,2)=0.5;
%mapC(:,:,3)=0.5;
    
ss=surf(pltx,plty,pltz,mapC)
set(ss,'LineStyle','none');


wf_porb_neg=wf_porb;
wf_porb_neg(wf_porb_neg>0)=0;
wf_porb_neg=abs(wf_porb_neg).^2;

wf_porb_plt=wf_porb_neg;
pltx0=sin(mesh1).*cos(mesh2).*wf_porb_plt;
plty0=sin(mesh1).*sin(mesh2).*wf_porb_plt;
pltz0=cos(mesh1).*wf_porb_plt;


pltx=pltx0*newx(1)+plty0*newy(1)+pltz0*newz(1)+shift_center(1);
plty=pltx0*newx(2)+plty0*newy(2)+pltz0*newz(2)+shift_center(2);
pltz=pltx0*newx(3)+plty0*newy(3)+pltz0*newz(3)+shift_center(3);



mapC=zeros(size(plty,1),size(plty,2),3);
mapC(:,:,3)=1;
    
ss=surf(pltx,plty,pltz,mapC)
set(ss,'LineStyle','none');


end